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Abstract 

The cross gramian matrix is a tool for model reduction and 
system identification, but it is only computable for square con¬ 
trol systems. For symmetric systems the cross gramian pos¬ 
sesses a useful relation to the system’s associated Hankel sin¬ 
gular values. Yet, many real-life models are neither square nor 
symmetric. In this work, concepts from decentralized control 
are used to approximate a cross gramian for non-symmetric and 
non-square systems. To illustrate this new non-symmetric cross 
gramian, it is applied in the context of model order reduction. 
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1 Introduction 

The cross gramian was introduced in for single-input-single-output (SISO) 
systems and extended to multiple-input-multiple-output (MIMO) systems in 
m- With many applications in model order reduction (MOR) and system 
theory in general, such as system identification m, decentralized control m 
parameter identification m or sensitivity analysis |20j . a major hindrance in 
the use of the cross gramian matrix is the constraint that it can be computed 
strictly for square systems and exhibits its core property only for symmetric 
systems. An extension of the cross gramian to non-symmetric systems enables 
such uses and particularly the (approximate) balancing state-space reduction by 
the computation of a single gramian matrix. This work can be seen as a follow 
up to m and |3j, which previously expanded the scope of cross gramian. 

The object of interest in this context is a linear time-invariant state-space 
system: 


x{t) = Ax{t) + Bu{t), 

y{t) = Cx{t) + Du{t), (1) 

x(0) = xo, 

^Contact: christian.hiinpe@uni-inuenster.de, inario.ohlberger@uni-niuenster.de, Insti¬ 
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which consists of a linear vector field composed of a system matrix A G 
and an input matrix B G as well as a linear output functional com¬ 
posed of an output matrix C G and a feed-forward matrix D G 


In the scope of this work, only a trivial feed-forward matrix 13 = 0 is consid¬ 
ered. Furthermore, the system is assumed to be asymptotically stable, hence all 
eigenvalues of the system matrix A lie in the left half-plane Re(Ai(A)) < 0. The 
classic cross gramian is defined for a subset of these systems and is expanded to 
arbitrary systems in the remainder of this work, which is organized as follows. 
In [Section'^ the cross gramian is reviewed. Existing approaches and the new 
result for the non-symmetric cross gramian are presented in |Section 3| Lastly, in 
[Section "dl verification and validation for the new non-symmetric cross gramian 
and comparison with established methods is conducted. 

2 The Cross Gramian 

With the controllability operator C{u) := e^* Bu{t) dt and the observability 

operator 0{x) := Ce^*' x, the controllability and observability of a system can 
be evaluated through the associated controllability gramian Wc ■= CC* and ob¬ 
servability gramian Wq '■= 0*0. A third system gramian, called cross gramian, 
combines controllability and observability into a single matrix. 

2.1 Square Systems 

For a square system, a system with the same number of inputs and outputs, 
the cross gramian is defined as the product of controllability operator C and 
observability operator O 0: 



Classically, the cross gramian is computed as solution to the Sylvester equation: 

AWx + WxA = -BC, 


which relates to the definition in ([^ through integration by parts: 



AWx = —BC — WxA. 
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2.2 Symmetric Systems 

A system in the form of (fTl) is called symmetric if the system’s gain is symmet- 

ri(E] 


CA-^B = {CA-^Bf. 

In other words for a symmetric system exists a symmetric matrix P such that: 


AP = PA^ and B = p-^C'^ 


( 3 ) 


are fulfilled [HJE]. Since a SISO system has a scalar gain, all SISO systems are 
symmetric. 

For a symmetric system, the absolute value of the eigenvalues of the cross 
gramian are equal to the Hankel singular values [19] : 


= WcWo 


^ |A(VFx)| = VHWcWo). 


( 4 ) 


This core property of the cross gramian allows to evaluate controllability and 
observability information of a system by computing a single gramian matrix. 
Numerically, this can be more efficient than computing two, namely the con¬ 
trollability and observability gramians, which additionally have to be balanced 
to be used, in example, for model order reduction through balanced truncation 


m- 


3 The Non-Symmetric Cross Gramian 

In this section, existing approaches for cross gramians of non-symmetric systems 
and selected methods from decentralized control are briefly summarized; the 
latter is employed to expand the scope of the cross gramian from symmetric to 
non-symmetric systems. 

3.1 Previous Work 

To the authors’ best knowledge there exist two methods to broaden the scope of 
the cross gramian for non-symmetric MIMO systems and towards more general 
configurations. 

The first approach from [I] extends the applicability of the cross gramian 
from symmetric systems to the wider class of orthogonally symmetric systems. 
Given a symmetric matrix P — for which AP = PA^ holds and an orthog¬ 
onal matrix U, with the property B — PCU^ if O < M, or C = PBU"^ if 
M < O, then the system is orthogonally symmetric and the associated cross 
gramian: 



satisfies the core property Q. 

^Equivalently the symmetry of the impulse response, transfer function or Markov parameter 
can be used. 
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The second approach, presented in [T5J[T3], uses embedding of a non-symmetric 
or non-square system into a symmetric system, and relies on a symmetrizer ma¬ 
trix [3] -/ = 

AJ = JA^. 


For a symmetrizer matrix J to A, an embedding is given by: 

x{t) = Ax{t) + i^JC^ B) u{t), 


y{t) = (^T^-i) 


the associated cross gramian has the form: 


Wx= I e^*{JC'^C + BB^J-^)e^*dt. 


T 


(5) 


While the first approach preserves the central property from Q in Wx for 
orthogonally symmetric systems, the second approach approximates it in Wx 
for arbitrary systems by an embedding. 


3.2 System Gramian Decomposition 

The basis for the non-symmetric cross gramian is a result from decentralized 
control, which aims to partition MIMO systems into sets of SISO systems with 
input-output pairings that exhibit the strongest coherence. To this end a rela¬ 
tion between the MIMO system gramians and the associated SISO subsystem 
gramians is described. As a first step a MIMO system is decomposed into 
M X O SISO systems, by partitioning the input matrix B and output matrix C 
column-wise and row-wise respectively: 


B = {hi ... bx!) ,bie 


nAfxl 


, C = 


Cl 


AO/ 


,Cj G 


nix Af 


Each combination of hi and cj induces a SISO system with the following system 
gramians: 


wh-.= 


pOO 

.— / 6 6 Clc, 

^0 


wl>-.= 


e'^ * cJ Cj e"^* dt. 




3^* biCj e'^* dt. 


The system gramians computed for the SISO subsystems relate to the full 
MIMO gramians as shown in m- 


M 

O 

Wo = Y/^o- 


( 6 ) 
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For square systems, also the cross gramian can be computed and the following 
identity holds m- 

M(=0) 

Wx= 

i=l 


3.3 Main Result 


Next, the previous result is utilized for the computation of an approximate 
cross gramian for non-square or non-symmetric systems. The central idea is to 
exploit, that for any SISO system a cross gramian with the property Q can be 
computed. 

With the product of controllability and observability can be expressed 
as sum of squared cross gramians: 


MO MO 

WcWo = 

2=1 j = l i—1 j — 1 


( 8 ) 


Due to the squaring of this ansatz is not numerically efficient. There¬ 

fore, an alternative approximate non-symmetric cross gramian, related to 0. 
is introduced. 

Definition 1 (Non-Symmetric Cross Gramian) 

The non-symmetric cross gramian is defined as the sum of the cross gramians 
of all M X O SISO subsystems: 


M O 

i=i i=l 


(9) 


Obviously, this gramian does not preserve the cross gramian’s property Q 
and it should be emphasized that the non-symmetric cross gramian does not 
reduce to the classic cross gramian in case of a symmetric system (compare 
0). Yet, for a linear system, Wz yields the following representation: 


M O 

i=lj=l■ 


3"^* biCj e^* dt 


, M O 
e"^* biCj dt 

i=i j=i 

, M O 

dt. 

i=l 3=1 


( 10 ) 


Hence, this approximate cross gramian Wz is equal to the cross gramian of the 
SISO system given by A, the row sum of the input matrix B and the column 
sum of the output matrix C, and thus can be seen as an “average” cross gramian 
over all SISO subsystems. 

Both approaches in [Section 3.1| share the common drawback, that they are 
limited to linear systems 0, and additionally may require a, potentially com¬ 
putationally expensive, symmetrizer. The new approach, proposed in 0 , does 
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not require the linear structure of the underlying system or derive system prop¬ 
erties using linear algebra. Only a cross gramian of a SISO system needs to 
be computable. Thus, the non-symmetric cross gramian can even be computed 
for nonlinear systems if a nonlinear cross gramian is available; this could be for 
example an empirical cross gramian m- Empirical gramians m are computed 
from (simulated) trajectories of the underlying system with perturbations in in¬ 
put and initial state. The empirical cross gramian has been introduced for SISO 
systems in EO] and extended to MIMO systems in m-, and as shown in m, the 
empirical cross gramian is equal to the cross gramian (|^ for linear systems Q , 
hence the empirical cross gramian can be used for the following experiments. 

4 Numerical Results 

The presented method is implemented as part of the empirical gramian frame¬ 
work (emgr) [5] introduced in [U]; and the following numerical experiments are 
conducted using emgr, which is compatible with OCTAVE and MATLAB®. 
The source code for reproducing the experiments can be found under an open 
source license in the supplemental materials and at http://runmycode.org/ 
companion/view/913 . 

From a computational point of view, the proposed method has the advan¬ 
tage, that for all SISO subsystem cross gramians the trajectories for perturbed 
initial states (observability), which consume the dominant fraction of overall 
computational time, have to be computed only once. 

Next, the presented non-symmetric cross gramian is tested in the context of 
projection-based model reduction [5]: For a linear system (Q a reduced order 
model is obtained through projection-based model reduction using the (trun¬ 
cated) projection Ui,Vi by: 

Xr = ViAUiXr{t) + ViBu{t), 

Ur = CU\Xr{t). 

Among others, such projections can be computed by balancing approaches, like 
balanced truncation m utilizing the controllability gramian and observabil¬ 
ity gramian, or direct truncation (approximate balancing) utilizing the cross 
gramian. By truncating the left singular vectors from the singular value decom¬ 
position of the cross gramian a (one-sided) Galerkin projection is generated, for 
further details see m, 

Wx UDV ^U={Ui U2) Vi := . 

Generally, a Galerkin projection cannot be expected to be stability-preserving, 
yet in case of the cross gramian the same argument as for balanced POD in EH 
Gh. 5.4] can be made due to the use of an energy-preserving inner-product, 

WxW^=COiCOy = (C,C*)oo-- 


This property transfers also to the newly introduced non-symmetric cross gramian, 
since it reduces to a classic cross gramian for the “averaged” system (101 that 
has the same system matrix A. 
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Figure 1: Relative L 2 output error of reduced order models from balanced 
truncation, cross gramian and non-symmetric cross gramian for a state-space 
symmetric system. 


First, the approximate non-symmetric cross gramian ® is compared to 
balanced truncation unm and the classic cross gramian (S for a symmetric 
system. 

A state-space symmetric system A = , B = of state-space dimension 

N = dim(a:(t)) = 64 and input dimension M = dim(u(t)) = dim{y{t)) = 0 = 8 
is selected, with a negative Lehmer matrix as system matrix A, Aij := — p j) 
and uniformly random generated input matrix B = C'^. 


To confirm (10), the error between the proposed non-symmetric cross gramian 
Wz and the cross gramian Wx of the SISO system {A, b, c}, with bi = J2j=i 
and Cj = i compared in the Frobenius norm: 

\\Wz-Wx\\f « 1.02E-13, 

which is close to machine precision in double precision floating point arithmetic. 

For the following experiments a zero initial state ccq = 0 is set and an impulse 
input = 5(t), i = 1... M, is applied. The relative L 2 output error is then 
evaluated for varying reduced order state-space dimensions. Such an coherence 
assessment for a reduced order impulse response is meaningful since a linear 
system response is a convolution of the impulse response with an input function. 
Figure l] shows, that the reduced order models obtained by balanced trun¬ 


cation and the cross gramian exhibits the same behavior due to state-space 
symmetric nature of the system. The newly proposed non-symmetric cross 
gramian from |Definition 1| does not achieve the same accuracy, but provides a 
lower output error for lesser order reduced models. A lower accuracy is to be 
expected due to the use of a one-sided projection; while balanced truncation 
uses a two-sided Petrov-Galerkin projection, which in the state-space symmet¬ 
ric case is equal to the projection obtained from the (symmetric) classic cross 
gramian. Notably, the output error of the reduced order model from the non- 
symmetric cross gramian drops already at a very low order n > 8 to the steady 
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Figure 2: Relative L 2 output error of reduced order models from balanced 
truncation, embedding cross gramian and non-symmetric cross gramian for a 
non-square system. 

error level while balanced truncation and the classic cross gramian reach this 
error not until a reduced order of n > 38. 

Second, the non-symmetric cross gramian is compared to balanced trunca¬ 
tion and the approximate cross gramian (|^ of the symmetric system derived by 
embedding from (§ for a non-square (and thus non-symmetric) system. 

To prevent the computation of a symmetrizer matrix J, the symmetric sys¬ 
tem matrix A S ]^64x64 example is reused, yet now, a uniformly 

random generated input matrix B G and a uniformly random output ma¬ 
trix C € is selected, thus the system is non-square and non-symmetric, 

since M = 4 and 0 = 8. Also for this example, zero initial state Xq = 0 and 
impulse input Ui{t) = S{t), i = 1...4, is applied and the relative L 2 output 
error is evaluated for varying reduced order state-space dimensions in [Figure 2| 

With the reference of balanced truncation, the cross gramian of the embed¬ 
ded system performs worse as predicted in dHlIIH]. Again, the non-symmetric 
cross gramian yields reduced models with less relative output error for small 
(reduced) orders. 

Lastly, a square but non-symmetric stable system with system matrix A G 
R64x 64^ uniformly random input matrix B G R^^xs uniformly random out¬ 
put matrix C G jg chosen. For this system, balanced truncation, the 

cross gramian and the non-symmetric cross gramian are compared. Due to the 
non-symmetric nature of the system the cross gramian has no theoretical foun¬ 
dation, yet is applicable since the system is square. Heuristically, also in [T] 
workable results for non-symmetric systems have been obtained for the classic 
cross gramian. 

[Figure 3| shows similar results as for the non-square system: the cross gramian 
and non-symmetric cross gramian both a lower accuracy, yet the non-symmetric 
cross gramian requires the least states for this error level even compared to bal¬ 
anced truncation. 

















Reduced State Space Dimension 

Figure 3: Relative L 2 output error of reduced order models from balanced 
truncation, cross gramian and non-symmetric cross gramian for a non-symmetric 
system. 

For all three experiments, the state-space model reduction error of the pre¬ 
sented non-symmetric cross gramian is worse than for balanced truncation, but 
the descent of the error is steeper. Hence, if the lower accuracy is acceptable, 
a smaller reduced order model can be constructed with this (non-symmetric) 
variant of the cross gramian method. This superior performance of the non- 
symmetric cross gramian compared to balanced truncation for low order reduced 
order models requires further investigation. 

It should be noted, that the non-symmetric cross gramian is not suitable 
for frequency-space-based approximation. This is due to its construction by 
averaging the sub-system cross gramians. For example. Hardy 7^2- and Hoo- 
errors between full and reduced order model will decay very slowly or not at all, 
hence this method is targeted purely at state-space approximations. 

5 Conclusion 

In this work a non-symmetric cross gramian, based on concepts from decen¬ 
tralized control, is proposed, and demonstrated to provide viable results for 
linear non-symmetric and non-square systems, which are outside the scope of 
the regular cross gramian. 

Future work will evaluate the effectiveness of the non-symmetric cross gramian 
for parametric and nonlinear systems. Furthermore, since the procedure sug¬ 
gested in m does not yield stable reduced order models in this setting, an 
investigation of alternative two-sided projections [7] for the (non-symmetric) 
cross gramian, may yield errors comparable to balanced truncation. 


9 















Acknowledgement 

This work was supported by the Deutsche Forschungsgemeinschaft, the Open 
Access Publication Fund of the University of Munster, DFG EXC 1003 Cells in 
Motion - Cluster of Excellence, Munster, Cermany as well as by the Center for 
Developing Mathematics in Interaction, DEMAIN, Munster, Cermany. 


References 

[1] U. Baur and P. Benner. Cross-gramian based model reduction for data- 
sparse systems. Electronic Transactions on Numerical Analysis^ 31:256- 
270, 2008. URL: http://eudinl.org/doc/130548 

[2] U. Baur, P. Benner, and L. Feng. Model order reduction for linear and non¬ 
linear systems: a system-theoretic perspective. Archives of Computational 
Methods in Engineering^ 21(4), 2014. doi; 10.1007/sll831-014-9111-2 

[3] K. Datta. The matrix equation XA — BX = R and its applications. 
Linear Algebra and its Applications, 109:91-105, 1988. doi: 10.1016/ 
0024-3795(88)90200-5 

[4] J.A. De Abreu-Carcia and F. Fairman. A note on cross grammians for or¬ 
thogonally symmetric realizations. IEEE Transactions on Automatic Con¬ 
trol, 31(9):866-868, 1986. doi:10.1109/TAC.1986.1104421 

[5] K.V. Fernando and H. Nicholson. On the structure of balanced and other 
principal representations of siso systems. IEEE Transactions on Automatic 
Control, 28(2):228-231, 1983. doi: 10.1109/TAC. 1983.1103195 

[6] L. Fortuna, A. Callo, C. Cuglielmino, and C. Nunnari. On the solution of 
a nonlinear matrix equation for mimo symmetric realizations. Systems & 
control letters, ll(l):79-82, 1988. doi: 10.1016/0167-6911(88)90115-6 

[7] S. Gugercin and A.C. Antoulas. An 7/2 error expression for the lanczos 
procedure. In ./.Snd IEEE Conference on Decision and Control, volume 2, 
pages 1869-1872, 2003. doi:10.1109/CDC.2003.1272886 

[8] C. Himpe. emgr - Empirical Gramian Framework (Version 3.5). http: 
//gramian.de, 2015. 

[9] C. Himpe and M. Ohlberger. A unified software framework for empirical 
gramians. Journal of Mathematics, 2013:1-6, 2013. doi: 10.1155/2013/ 
365909, 

[10] C. Himpe and M. Ohlberger. Cross-gramian based combined state and 
parameter reduction for large-scale control systems. Mathematical Problems 
in Engineering, 2014:1-13, 2014. doi: 10.1155/2014/843869, 

[11] P. Holmes, J.L. Lumley, G. Berkooz, and C.W. Rowley. Turbulence, Co¬ 
herent Structures, Dynamical Systems and Symmetry. Cambridge Mono¬ 
graphs on Mechanics. Cambridge University Press, 2012. doi: 10.1017/ 
CB09780511919701, 


10 








[12] S. Lall, J.E. Marsden, and S. Glavaski. Empirical model reduction of con¬ 
trolled nonlinear systems. Proceedings of the IFAC World Congress, F:473- 
478,1999. URL:http://authors.library.caltech.edu/20343/2/10.1. 
1.123.4669.pdf 

[13] A.J. Laub, L.M. Silverman, and M. Verma. A note on cross-grammians 
for symmetric realizations. Proceedings of the IEEE, 71(7):904-905, 1983. 
doi:10.1109/PR0C.1983.12688 

[14] L.A. Mironovskii and T.N. Solv’eva. Analysis of multiplicity of han- 
kel singular values of control systems. Automation and Remote Control, 
76(2):205-218, 2015. doi:10.1134/S0005117915020022, 

[15] B. Moaveni and A. Khaki-Sedigh. Input-output pairing based on cross- 
gramian matrix. International Joint Conference SICE-ICAS, pages 2378- 
2380, 2006. doi:10.1109/SICE.2006.314989 

[16] B. Moore. Principal component analysis in linear systems: Controllabil¬ 
ity, observability, and model reduction. IEEE Transactions on Automatic 
Control, 26(l):17-32, 1981. doi: 10.1109/TAC. 1981.1102568 

[17] S. Rahrovani, M.K. Vakilzadeh, and T. Abrahamsson. On gramian-based 
techniques for minimal realization of large-scale mechanical systems. In 
Topics in Modal Analysis, Volume 1, pages 797-805. Springer, 2014. doi: 
10.1007/978-l-4614-6585-0_75 

[18] D.C. Sorensen and A.C. Antoulas. Projection methods for balanced model 
reduction. Technical Report, pages 1-18, 2001. URL: http://www. caam. 
rice . edu/caain/trs/2001/TR01-03. pdf 

[19] D.C. Sorensen and A.C. Antoulas. The Sylvester equation and approximate 
balanced reduction. Linear Algebra and its Applications, 351:671-700, 2002. 
doi:10.1016/30024-3795(02)00283-5 

[20] S. Streif, R. Findeisen, and E. Bullinger. Relating cross gramians 
and sensitivity analysis in systems biology. Theory of Networks and 
Systems, 10.4:437-442, 2006. URL: http://eprints.nuim. ie/1768/1/ 
HamiltonGramian.pdf 

[21] J.C. Willems. Realization of systems with internal passivity and symmetry 
constraints. Journal of the Eranklin Institute, 301(6):605-621, 1976. doi: 
10.1016/0016-0032(76)90081-8 


11 



